################################################################################
class rap_rad(rap_trans):
  """ 
  Class for storing Einstein A coefficients.
  """
  def __init__(self,filename=None,levs=None,**args):
    self.clear()
    if levs is not None: self.set_levs(levs)
    if filename is not None: self.loadfile(filename,**args)

# ------------------------------------------------------------------------------

  def clear(self,keeplevs=False):
    """ Erase all data. """
    rap_trans.clear(self,keeplevs=keeplevs)
    self.v_lifetime={}

# ------------------------------------------------------------------------------

  def loadfile(self,filename,cup=None,iminlev=None,imaxlev=None,\
                                      fminlev=None,fmaxlev=None):
    """ Load radiative data from file. """
    ftype=rlf_filetype(filename)
    if ftype == 'olg':
      self.read_olg(filename,cup=cup,iminlev=iminlev,imaxlev=imaxlev,\
                                     fminlev=fminlev,fmaxlev=fmaxlev)
    elif ftype == 'adf04':
      self.read_adf04(filename,cup=cup,iminlev=iminlev,imaxlev=imaxlev,\
                                       fminlev=fminlev,fmaxlev=fmaxlev)
    elif ftype == 'chianti_wgfa':
      self.read_chianti(filename,iminlev=iminlev,imaxlev=imaxlev,\
                                 fminlev=fminlev,fmaxlev=fmaxlev)
    elif ftype == 'x050':
      self.read_xstar(filename,iminlev=iminlev,imaxlev=imaxlev,\
                                fminlev=fminlev,fmaxlev=fmaxlev)
    elif ftype == 'nist_lines':
      self.read_nist(filename,iminlev=iminlev,imaxlev=imaxlev,\
                              fminlev=fminlev,fmaxlev=fmaxlev)
    elif ftype == 'tiptop_lines':
      self.read_tiptop(filename,iminlev=iminlev,imaxlev=imaxlev,\
                                fminlev=fminlev,fmaxlev=fmaxlev)
    else:
      raise ErrorRAP('unsupported file type')

# ------------------------------------------------------------------------------

  def reindex(self,remap):
    """ 
    Reassign level indices.  Indices are remapped according to dictionary,
    remap, old (key) -> new (value).  
    """
    imap=rap_trans.reindex(self,remap)
    self.v_lifetime=dict([(imap[i],val) for (i,val) in self.v_lifetime.items()])
    return imap

# ------------------------------------------------------------------------------

  def avalue(self,i,j=None,imap=False):
    """ 
    Return the A-value for the i-j transition.
    
    Note: it is assumed that state i has a lower energy than state j.
    """
    return self.data(i,j,imap)

# ------------------------------------------------------------------------------

  def set_avalue(self,i,j,avalue,imap=False):
    """ Add or update an A-value (also updates lifetime). """
    ii,jj=i,j
    if imap: (ii,jj)=self.v_levs.lmap((ii,jj))
    if self.v_levs.energy(ii) > self.v_levs.energy(jj): ii,jj=jj,ii
    (iimin,iimax)=self.initial_limit()
    if ii < iimin and ii > iimax:
      msg='initial index outside of user-defined range: %s'%ii
      raise ErrorRAP(msg)
    (jjmin,jjmax)=self.final_limit()
    if jj < jjmin and jj > jjmax:
      msg='final index outside of user-defined range: %s'%jj
      raise ErrorRAP(msg)
    if (ii,jj) in self:
      tau=self.v_lifetime[jj]
      taup=1./(1./tau-self[(ii,jj)]+avalue)
    else:
      taup=1./avalue
    self.v_data[(ii,jj)]=avalue
    self.v_lifetime[jj]=taup

# ------------------------------------------------------------------------------

  def set_avalues(self,avalues,imap=False):
    """ 
    Efficiently add many new A-values.  Input is a dictionary with keys being
    the transition tuple by index (lower,upper) and the values are the A-values.
    """
    newsum={}
    oldsum={}
    for (i,j),av in avalues.iteritems():
      ii,jj=i,j
      if imap: (ii,jj)=self.v_levs.lmap((ii,jj))
      if self.v_levs.energy(ii) > self.v_levs.energy(jj): ii,jj=jj,ii
      (iimin,iimax)=self.initial_limit()
      if ii < iimin and ii > iimax:
        msg='initial index outside of user-defined range: %s'%ii
        raise ErrorRAP(msg)
      (jjmin,jjmax)=self.final_limit()
      if jj < jjmin and jj > jjmax:
        msg='final index outside of user-defined range: %s'%jj
        raise ErrorRAP(msg)
      if jj not in newsum: 
        newsum[jj]=0.0
        oldsum[jj]=0.0
      newsum[jj]+=av
      if (ii,jj) in self: oldsum[jj]+=self[(ii,jj)]
      self.v_data[(ii,jj)]=av
    for jj in newsum:
      if jj not in self.v_lifetime:
        self.v_lifetime[jj]=1./newsum[jj]
      else:
        self.v_lifetime[jj]=1./(1./self.v_lifetime[jj]-oldsum[jj]+newsum[jj])


# ------------------------------------------------------------------------------

  def gavalue(self,i,j=None,imap=False):
    """ 
    Return the gA-value for the i-j transition.
    
    Note: it is assumed that state i has a lower energy than state j.
    """
    ii,jj=i,j
    if jj is None and isinstance(ii,tuple): (ii,jj)=ii
    a=self.avalue(ii,jj,imap=imap)
    g=self.levs().weight(jj,imap=imap)
    return a*g

# ------------------------------------------------------------------------------

  def gfvalue(self,i,j=None,imap=False):
    """ 
    Return the gf-value for the i-j transition.
    
    Note: it is assumed that state i has a lower energy than state j.
    """
    ii,jj=i,j
    if jj is None and isinstance(ii,tuple): (ii,jj)=ii
    ga=self.gavalue(ii,jj,imap=imap)
    factor=abs(self.levs().energy(ii,imap=imap)-\
               self.levs().energy(jj,imap=imap))
    factor=0.66702*(factor*INVERSECM)**2
    return ga/factor

# ------------------------------------------------------------------------------

  def lifetime(self,i,imap=False):
    """ Return the lifetime of level, i.  Return None if infinite. """
    ii=i
    if imap:
      ii=self.v_levs.lmap(ii)
      if ii is None: raise ErrorRAP('invalid mapping')
    if ii not in self.v_lifetime: return None
    return self.v_lifetime[ii]

# ------------------------------------------------------------------------------

  def prune(self,fraction=0.001):
    """ 
    Remove A-values which compose less than a fraction (default: 0.001) of 
    the (inverse) lifetime of the upper state.
    """
    trans=self.v_data.keys()
    for (i,j) in trans:
      chk=1./self.v_lifetime[j]
      if self.v_data[(i,j)] < fraction*chk: del self.v_data[(i,j)]

# ------------------------------------------------------------------------------

  def set_data(self,i,j=None,value=None,imap=False):
    """ 
    Set data value for the i-j transition.  This is a replacement for the
    set_data() function in rap_trans, as to allow for the lifetime to be 
    properly adjusted.
    
    Note: it is assumed that state i has a lower energy than state j.
    """
    if isinstance(i,tuple):
      (ii,jj)=i
      tval=j
    else:
      ii,jj=i,j
      tval=value
    if tval is None:
      raise ErrorRAP("need to give a new value")
    if ii not in self.levs():
      raise ErrorRAP("Invalid state index given: %s"%ii)
    if jj not in self.levs():
      raise ErrorRAP("Invalid state index given: %s"%jj)
    trans=(ii,jj)
    if imap:
      trans=self.v_levs.lmap(trans)
      if trans is None: raise ErrorRAP('invalid mapping')

    oldval=0.0
    if trans in self: oldval=self[trans]
    self.v_data[trans]=tval

    # need to adjust lifetime
    (ii,jj)=trans
    tau=self.lifetime(jj)
    if tau is None:
      taup=1./tval
    else:
      taup=1./(1./tau-oldval+tval)
    self.v_lifetime[jj]=taup

# ------------------------------------------------------------------------------

  def LS_average(self,omit_incomplete=True):
    """
    Recast LSJ radiative data into LS coupling.  New values are averaged over
    the initial (upper) states and summed over final (lower) states.  The
    option, omit_incomplete, allows the skipping of terms which do not have all
    fine-structure levels present. 
    """
    if not isinstance(self.v_levs,rap_states_lsj):
      raise ErrorRAP("radiative data need to be in LSJ coupling")

    levs_ls=self.v_levs.LS_average(omit_incomplete=omit_incomplete)


    # for each LSJ state, get weighting factors for the constituent levels
    tmap={}        # level -> term map
    wmap={}        # weighting factor for each level
    for idx in levs_ls:
      (conf,term,alfa,par,e)=levs_ls.level(idx)
      wsum=float(rlf_termweight(term))

      tidxs=self.v_levs.search(conf=conf,term=term,alpha=alfa)
      for tidx in tidxs:
        wgt=self.v_levs.weight(tidx)
        wmap[tidx]=float(wgt)/wsum
        tmap[tidx]=idx

    # average rates
    toadd={}
    for (i,j),aval in self.v_data.iteritems():
      t=(tmap[i],tmap[j])
      if t not in toadd: toadd[t]=0.0
      toadd[t]+=wmap[j]*aval

    # new object
    rad=rap_rad(levs=levs_ls)
    rad.v_lifetime=dict([(i,0.0) for i in rad.v_levs.indices()])
    for (i,j),aval in toadd.iteritems():
      rad.v_data[(i,j)]=aval
      rad.v_lifetime[j]+=aval

    for j in rad.v_lifetime.keys():
      asum=rad.v_lifetime[j]
      if asum == 0.0: 
        del rad.v_lifetime[j]
      else:
        rad.v_lifetime[j]=1.0/asum

    return rad

# ------------------------------------------------------------------------------

  def LSJ_split(self,levs_lsj=None):
    """ 
    Split LS energy radiative data in LSJ states.  Need to pass in the LSJ
    level structure (ideally calculated using rap_states_ls.LSJ_split() )
    which will be the level structure associated with the new rap_rad object.
    Dipole allowed transitions will be split according to Eq. 14.50 of Cowan
    (1981).  This assumes that LS coupling is a good description of the system
    and the dipole contribution to the A-value is dominant.  Non-dipole
    transitions will be split statistically and cannot be trusted.  Ideally,
    this routine will serve as a place holder for more accurate values to be
    read in.
    """

    # take term index and get list of level indices (levs_lsj)
    def mapme(i):
      conf=self.levs().configuration(i)
      term=self.levs().term(i)
      alpha=self.levs().seniority(i)
      tlist=levs_lsj.cfilter(conf)
      tlist=levs_lsj.tfilter(term,alpha=alpha,subset=tlist)
      return tlist

    if not isinstance(self.v_levs,rap_states_ls):
      raise ErrorRAP("radiative data need to be in LS coupling")
    if not isinstance(levs_lsj,rap_states_lsj):
      raise ErrorRAP("need to supply rap_states_lsj instance")

    toadd={}
    lmap={}
    for (ilow,iupp) in self:
      if ilow not in lmap: lmap[ilow]=mapme(ilow)
      if iupp not in lmap: lmap[iupp]=mapme(iupp)

      # check if dipole allowed
      dip=True
      (glow,llow)=rlf_termdecode(self.levs().term(ilow))
      (gupp,lupp)=rlf_termdecode(self.levs().term(iupp))
      if glow != gupp: dip=False
      if (llow-lupp) not in [-1,0,+1]: dip=False
      if llow == 0 and lupp == 0: dip=False

      slow=0.5*float(glow-1)
      supp=0.5*float(gupp-1)

      wgts={}
      for ifs_low in lmap[ilow]:
        jlow=levs_lsj.j(ifs_low)

        for ifs_upp in lmap[iupp]:
          jupp=levs_lsj.j(ifs_upp)

          if dip:     # Cowan 14.50
            wgt=(2.*jlow+1.)*(2.*jupp+1.)*\
                pyang_wigner6j(llow,slow,jlow,jupp,1.,lupp)**2
          else:
            wgt=(2.*jlow+1)
          wgts[(ifs_low,ifs_upp)]=wgt

      twgt=sum(wgts.values())
      factor=1./twgt
      aval=self.avalue(ilow,iupp)
      for key in wgts: toadd[key]=factor*wgts[key]*aval

    # assemble new object
    out=rap_rad()
    out.set_levs(levs_lsj)
    out.v_lifetime=dict([(i,0.0) for i in levs_lsj.indices()])
    for (ilow,iupp),aval in toadd.iteritems(): 
      out.v_data[(ilow,iupp)]=aval
      out.v_lifetime[iupp]+=aval
    for j in out.v_lifetime.keys():
      asum=out.v_lifetime[j]
      if asum == 0.0: 
        del out.v_lifetime[j]
      else:
        out.v_lifetime[j]=1.0/asum
    return out

# ------------------------------------------------------------------------------

  def read_olg(self,filename,cup=None,iminlev=None,imaxlev=None,
                                      fminlev=None,fmaxlev=None):
    """ Read radiative data from olg file """
    if self.levs() == None: self.v_levs=rlf_loadstates(filename,cup=cup)
    tcup=self.v_levs.coupling()
    if tcup == 'LSJ':
      count=2
      j0=9
      j1=18
      i0=18
      i1=22
      ae0=33
      ae1=42
      am0=48
      am1=57
    elif tcup == 'LS':
      count=1
      j0=5
      j1=13
      i0=13
      i1=17
      ae0=23
      ae1=32
    else:
      msg='coupling scheme must be LS or LSJ for olg files'
      raise ErrorRAP(msg)

    self.clear(keeplevs=True)
    self.set_initial_limit(iminlev,imaxlev)
    self.set_final_limit(fminlev,fmaxlev)
    self.v_lifetime=dict([(i,0.0) for i in self.v_levs.indices()])
    processing=False
    fin=open(filename,'r')
    for line in fin:
      if 'E0-DATA' in line: count-=1
      if count < 0: break
      if count > 0: continue

      if not processing:
        if '-DATA' in line:
          processing=True
          tmp=re.findall(r'E(\d+).*DATA',line)
          rank=int(tmp[0])
        continue

      if line.strip() == '':
        processing=False
        continue
      j=int(line[j0:j1])
      i=int(line[i0:i1])
      if fminlev and j < fminlev: continue
      if fmaxlev and j > fmaxlev: 
        processing=False
        continue
      try:
        a=float(line[ae0:ae1])     # fails when < 1E-99 (e.g. 6.2456-101)
      except:
        a=0.0
      if tcup == 'LSJ' and rank > 1: a+=float(line[am0:am1])
      self.v_lifetime[j]+=a
      if iminlev and i < iminlev: continue
      if imaxlev and i > imaxlev: continue
      trans=(i,j)
      if trans not in self.v_data: self.v_data[trans]=0.0
      self.v_data[trans]+=a
    fin.close()

    for j in self.v_lifetime.keys():
      asum=self.v_lifetime[j]
      if asum == 0.0: 
        del self.v_lifetime[j]
      else:
        self.v_lifetime[j]=1.0/asum

# ------------------------------------------------------------------------------

  def read_adf04(self,filename,cup=None,iminlev=None,imaxlev=None,
                                        fminlev=None,fmaxlev=None):
    """ Read A-values from adf04 file """
    if self.levs() == None: self.v_levs=rlf_loadstates(filename,cup=cup)
    self.clear(keeplevs=True)
    tcup=self.v_levs.coupling()
    if cup is not None and cup != tcup:
      msg='file coupling scheme, %s, incompatible with given: %s'%(tcup,cup)
      raise ErrorRAP(msg)
    self.set_initial_limit(iminlev,imaxlev)
    self.set_final_limit(fminlev,fmaxlev)
    self.v_lifetime=dict([(i,0.0) for i in self.v_levs.indices()])
    fin=open(filename,'r')
    while 1:
      line=fin.readline()
      if line.split()[0] == '-1': break
    line=fin.readline()
    line=fin.readline()
    avalues=[]
    tmp=line.split()
    j=tmp[0]
    i=tmp[1]
    a=tmp[2]
    icol1=line.index(j)+len(j)
    icol2=line[icol1:].index(i)+len(i)+icol1
    icol3=line[icol2:].index(a)+len(a)+icol2
    j=int(j)
    i=int(i)
    a=rlf_adas2float(a)
    if fminlev and j >= fminlev and fmaxlev and j <= fmaxlev:
      self.v_lifetime[j]+=a
      if iminlev and i >= iminlev and imaxlev and i <= imaxlev:
        trans=(i,j)
        self.v_data[trans]=a
    while 1:
      line=fin.readline()
      if line.strip() == '-1': break
      try:
        j=int(line[:icol1])
        i=int(line[icol1:icol2])
      except:
        break
      if (fminlev and j < fminlev) or (fmaxlev and j > fmaxlev): continue
      a=rlf_adas2float(line[icol2:icol3])
      self.v_lifetime[j]+=a
      if (iminlev and i < iminlev) or (imaxlev and i > imaxlev): continue
      trans=(i,j)
      self.v_data[trans]=a
    fin.close()

    for j in self.v_lifetime.keys():
      asum=self.v_lifetime[j]
      if asum == 0.0: 
        del self.v_lifetime[j]
      else:
        self.v_lifetime[j]=1.0/asum

# ------------------------------------------------------------------------------

  def read_chianti(self,filename,iminlev=None,imaxlev=None,
                                 fminlev=None,fmaxlev=None):
    """ Read A-values from CHIANTI wfga file """
    id_idx=0        # level index of lower state
    id_jdx=1        # level index of upper state
    id_wl=2         # wavelength (Angstroms)
    id_gf=3         # gf-value (1/s); can be zero if not given
    id_aval=4       # A-value (1/s)
    id_iconf=5      # lower state configuration
    id_istate=6     # lower state term and J-value in spectroscopic notation
    id_dash=7       # separator between upper and lower state information
    id_jconf=8      # upper state configuration
    id_jstate=9     # upper state term and J-value in spectroscopic notation

    if self.levs() == None: raise ErrorRAP('must provide levels object')
    self.clear(keeplevs=True)
    ne=self.levs().num_electrons()
    self.set_initial_limit(iminlev,imaxlev)
    self.set_final_limit(fminlev,fmaxlev)
    self.v_lifetime=dict([(i,0.0) for i in self.v_levs.indices()])

    indices=self.levs().indices()
    fin=open(filename,'r')
    while 1:
      line=fin.readline()
      tmp=line.split()
      if line.strip() == '' or int(tmp[0]) == -1: break
      if len(tmp) != 10: continue
      idx=int(tmp[id_idx])
      jdx=int(tmp[id_jdx])
      if imaxlev and idx > imaxlev: continue
      if fminlev and jdx < fminlev: continue
      if fmaxlev and jdx > fmaxlev: continue
      aval=float(tmp[id_aval])
      iconf=tmp[id_iconf]
      istate=tmp[id_istate]
      jconf=tmp[id_jconf]
      jstate=tmp[id_jstate]
      if idx not in indices or jdx not in indices:
        raise ErrorRAP('level indices not in level structure')

      cchk=self.levs().configuration(idx)
      tchk=self.levs().term(idx)+'%.1f'%(float(0.5*self.levs().j2(idx)))
      ic=rap_config(iconf,ne=ne)
      if ic != str(cchk) or istate != tchk:
        raise ErrorRAP('inconsistency with level labels')

      cchk=self.levs().configuration(jdx)
      tchk=self.levs().term(jdx)+'%.1f'%(float(0.5*self.levs().j2(jdx)))
      jc=rap_config(jconf,ne=ne)
      if jc != str(cchk) or jstate != tchk:
        raise ErrorRAP('inconsistency with level labels')

      self.v_lifetime[jdx]+=aval
      if iminlev and idx < iminlev: continue
      self.v_data[(idx,jdx)]=aval
    fin.close()

    for j in self.v_lifetime.keys():
      asum=self.v_lifetime[j]
      if asum == 0.0: 
        del self.v_lifetime[j]
      else:
        self.v_lifetime[j]=1.0/asum

# ------------------------------------------------------------------------------

  def read_tiptop(self,filename,iminlev=None,imaxlev=None,
                                fminlev=None,fmaxlev=None):
    """ 
    Read A-values from TIPTOP gf-value file.

    Note: It is assumed that the symmetry ordering in the TIPTOP files are
    ordered by energy.  That is, the iLV index increases with energy for each
    iSLP symmetry.  If this is not the case, then this routine needs to also
    take in the level file where the symmetries are explicitly defined and
    both the rap_states and rap_rad instances are created here.
    """
    fin=open(filename,'r')

    # get column info
    line=fin.readline()
    hline=fin.readline()
    line=fin.readline()
    tmp=hline.split()
    id_nz=tmp.index('NZ')
    id_ne=tmp.index('NE')
    id_islp=tmp.index('iSLP')
    id_jslp=tmp.index('jSLP')
    id_ilv=tmp.index('iLV')
    id_jlv=tmp.index('jLV')
    id_gf=tmp.index('gF')
    id_wl=tmp.index('WL(A)')

    if self.levs() == None: raise ErrorRAP('must provide levels object')
    self.clear(keeplevs=True)
    ne=self.levs().num_electrons()
    self.set_initial_limit(iminlev,imaxlev)
    self.set_final_limit(fminlev,fmaxlev)
    self.v_lifetime=dict([(i,0.0) for i in self.v_levs.indices()])

    # get mapping from symmetry to level index
    allterms=self.v_levs.terms()
    symmap={}
    for t in allterms:
      (g,l)=rlf_termdecode(t)

      # even
      sym="%s%s0"%(g,l)
      idxs=self.v_levs.symfilter(g,l,0)
      if len(idxs) > 0: symmap[sym]=idxs

      # odd
      sym="%s%s1"%(g,l)
      idxs=self.v_levs.symfilter(g,l,1)
      if len(idxs) > 0: symmap[sym]=idxs

    nz=None
    indices=self.v_levs.indices()
    while 1:
      line=fin.readline()
      tmp=line.split()
      if '=====' in line or line == '': break

      if nz is None:          # check that ion is compatible with levs
        nz=int(tmp[id_nz])
        ne=int(tmp[id_ne])
        if nz != self.v_levs.nuclear_charge() or \
           ne != self.v_levs.num_electrons():
          raise ErrorRAP("incompatible ion")

      islp=tmp[id_islp]
      ilv=int(tmp[id_ilv])
      jslp=tmp[id_jslp]
      jlv=int(tmp[id_jlv])
      gf=float(tmp[id_gf])

      # get indices, want idx to be lower level
      idx=symmap[islp][ilv-1]
      jdx=symmap[jslp][jlv-1]
      if idx not in indices or jdx not in indices:
        raise ErrorRAP('level indices not in level structure')
      if self.v_levs.energy(idx) > self.v_levs.energy(jdx): idx,jdx=jdx,idx
      if imaxlev and idx > imaxlev: continue
      if fminlev and jdx < fminlev: continue
      if fmaxlev and jdx > fmaxlev: continue

      g=self.v_levs.weight(jdx)
      s2=INVERSECM*(self.v_levs.energy(jdx)-self.v_levs.energy(idx))
      s2=s2**2
      aval=GF2A*s2*gf/g

      self.v_lifetime[jdx]+=aval
      if iminlev and idx < iminlev: continue
      self.v_data[(idx,jdx)]=aval
    fin.close()

    for j in self.v_lifetime.keys():
      asum=self.v_lifetime[j]
      if asum == 0.0: 
        del self.v_lifetime[j]
      else:
        self.v_lifetime[j]=1.0/asum

# ------------------------------------------------------------------------------

  def read_nist(self,filename,iminlev=None,imaxlev=None,
                              fminlev=None,fmaxlev=None):
    """ Read A-values from NIST text file. """
    id_obwl=0
    id_thwl=1
    id_aval=3
    id_ens=5
    id_confs=6
    id_terms=7
    id_js=8

    if self.levs() == None: raise ErrorRAP('must provide levels object')
    self.clear(keeplevs=True)
    ne=self.levs().num_electrons()
    self.set_initial_limit(iminlev,imaxlev)
    self.set_final_limit(fminlev,fmaxlev)
    self.v_lifetime=dict([(i,0.0) for i in self.v_levs.indices()])

    repchars=['+x','[',']','?']     # these can appear in NIST energy values

    fin=open(filename,'r')
    for line in fin:
      if '---' in line: continue
      if '?' in line: line=line.replace('?',' ')   # these can be annoying
      tmp=line.split('|')
      try:
        obwl=float(tmp[id_obwl])
      except:
        obwl=-1.0
      try:
        thwl=float(tmp[id_thwl])
      except:
        thwl=-1.0
      if obwl < 0.0 and thwl < 0.0: continue
      try:
        aval=float(tmp[id_aval])
      except:
        continue
      if tmp[id_ens].strip() == '': continue
      ttmp=tmp[id_ens].split('-')
      elow=ttmp[0]
      eupp=ttmp[1]
      for ch in repchars:
        elow=elow.replace(ch,'')
        eupp=eupp.replace(ch,'')
      elow=float(elow)/INVERSECM
      eupp=float(eupp)/INVERSECM
      ttmp=tmp[id_confs].split('-')
      clow=rap_config(ttmp[0],ne=ne)
      cupp=rap_config(ttmp[1],ne=ne)
      ttmp=tmp[id_terms].split('-')
      tlow=ttmp[0].replace('*','').strip()
      if '[' in tlow: tlow=None
      tupp=ttmp[1].replace('*','').strip()
      if '[' in tupp: tupp=None
      ttmp=tmp[id_js].split('-')
      tj2low=ttmp[0].strip()
      tj2upp=ttmp[1].strip()
      j2low=int(2*eval(tj2low+'.'))
      j2upp=int(2*eval(tj2upp+'.'))
      idx=self.levs().search(clow,tlow,j2low,energy=elow)
      jdx=self.levs().search(cupp,tupp,j2upp,energy=eupp)

      if fminlev and jdx < fminlev: continue
      if fmaxlev and jdx > fmaxlev: continue
      if imaxlev and idx > imaxlev: continue
      self.v_lifetime[jdx]+=aval
      if iminlev and idx < iminlev: continue
      self.v_data[(idx,jdx)]=aval
    fin.close()

    for j in self.v_lifetime.keys():
      asum=self.v_lifetime[j]
      if asum == 0.0: 
        del self.v_lifetime[j]
      else:
        self.v_lifetime[j]=1.0/asum

# ------------------------------------------------------------------------------

  def read_xstar(self,filename,iminlev=None,imaxlev=None,
                               fminlev=None,fmaxlev=None):
    """ Read A-values from xstar type 50 file """
    id_idx=9        # level index of lower state
    id_jdx=10       # level index of upper state
    id_wl=6         # wavelength (Angstroms)
    id_gf=7         # gf-value (1/s); can be zero if not given
    id_aval=8       # A-value (1/s)
    id_ionno=12     # xstar ion number

    if self.levs() == None: raise ErrorRAP('must provide levels object')
    self.clear(keeplevs=True)
    ne=self.levs().num_electrons()
    self.set_initial_limit(iminlev,imaxlev)
    self.set_final_limit(fminlev,fmaxlev)
    self.v_lifetime=dict([(i,0.0) for i in self.v_levs.indices()])

    fin=open(filename,'r')
    while 1:
      line=fin.readline()
      tmp=line.split()
      if line.strip() == '' or int(tmp[0]) == -1: break
      if tmp[-1] != '%' or int(tmp[0]) != 50:
        print line
        raise ErrorRAP('do not understand xstar type 50 file')
      idx=int(tmp[id_idx])
      jdx=int(tmp[id_jdx])
      if idx not in self.v_levs: continue
      if jdx not in self.v_levs: continue
      if fminlev and jdx < fminlev: continue
      if fmaxlev and jdx > fmaxlev: continue
      if imaxlev and idx > imaxlev: continue
      aval=float(tmp[id_aval])
      ionno=int(tmp[id_ionno])
      (z,ne)=xstar_iondecode(ionno)
      if z != self.levs().nuclear_charge() or ne != self.levs().num_electrons():
        raise ErrorRAP('ion in type 50 file does not match level structure')

      self.v_lifetime[jdx]+=aval
      if iminlev and idx > iminlev: continue
      self.v_data[(idx,jdx)]=aval
    fin.close()

    for j in self.v_lifetime.keys():
      asum=self.v_lifetime[j]
      if asum == 0.0: 
        del self.v_lifetime[j]
      else:
        self.v_lifetime[j]=1.0/asum

# ------------------------------------------------------------------------------

  def write_xstar(self,filename,subset=None,precision=5,imap=True):
    """ 
    Write xstar type 50 file. Must have ionization potential and nuclear charge
    specified in levels structure.
    """
    if self.levs().ionization_threshold() == None:
      raise ErrorRAP('must define ionization threshold')
    if self.levs().nuclear_charge() == None:
      raise ErrorRAP('must define nuclear charge')
    indices=self.levs().indices()
    if subset is not None: indices=subset
    trans=[]
    for (i,j) in self.transitions():
      if i in indices and j in indices: trans.append((i,j))

    zeff=self.levs().nuclear_charge()
    ne=self.levs().num_electrons()
    ionno=xstar_ionencode(zeff,ne)

    aform='%'+str(8+precision)+'.'+str(precision)+'E'
    iform='%6i'
    prefix=iform%50+iform%4+iform%0+iform%3+iform%4+iform%0
    suffix=iform%zeff+iform%ionno+' %'
    fout=open(filename,'w')
    for (i,j) in trans:
      try:
        wl=self.wavelength(i,j)
        gf=self.gfvalue(i,j)
      except:
        wl=0.0
        gf=0.0
      av=self.avalue(i,j)
      if av == 0.0: continue
      ii=i
      jj=j
      if imap:
        ii=self.v_levs.lmap_reverse(i)
        jj=self.v_levs.lmap_reverse(j)
        if ii is None or jj is None: continue
      line=prefix+aform%wl+aform%gf+aform%av+iform%ii+iform%jj+suffix
      fout.write(line+'\n')
    fout.close()


